Year of ALOS-PALSAR mosaic: 2019 Backscatter processing variant: maskLWUwb_med5_LCinterp AGB variant: cap999pctl_maskWUwb
# Histogram of external map ----
bacc_hist <- plot_agb_hist_density(agb_fps$bacc, agb_pal)
## [1] "Sampling..."
hist_fp <- here::here('data/reports/external_agb', str_glue('agb_hist_Baccini.png'))
dir.create(dirname(hist_fp))
# ggsave(hist_fp, width = 5, height = 3, dpi = 120)
(esa_hist <- plot_agb_hist_density(agb_fps$esa, agb_pal))
## [1] "Sampling..."
hist_fp <- here::here('data/reports/external_agb', str_glue('agb_hist_ESA.png'))
# ggsave(hist_fp, width = 5, height = 3, dpi = 120)
(glob_hist <- plot_agb_hist_density(agb_fps$glob, agb_pal, xlim=NULL))
## [1] "Sampling..."
hist_fp <- here::here('data/reports/external_agb', str_glue('agb_hist_GlobB.png'))
# ggsave(hist_fp, width = 5, height = 3, dpi = 120)
(avit_hist <- plot_agb_hist_density(agb_fps$avit, agb_pal, xlim=NULL))
hist_fp <- here::here('data/reports/external_agb', str_glue('agb_hist_Avit.png'))
# ggsave(hist_fp, width = 5, height = 3, dpi = 120)
# Load data
smmry_csv <- here::here(comparison_dir, '07_overview_by_LC.csv')
stats_by_lc <- read_csv(smmry_csv) %>%
filter(!is.na(name))
# Helpers
lc_tbl <- list(`4` = 'Tree cover', `6` = 'Shrubs', `5` = 'Grassland',
`3` = 'Bareland', `2` = 'Urban', `1` = 'Water')
(map_order <- stats_by_lc %>%
mutate(map_code = factor(map_code, levels = c('internal', 'esa', 'glob', 'avit', 'bacc'))) %>%
arrange(desc(map_code)) %>%
distinct(name) %>%
deframe())
## [1] "Baccini" "Avitabile"
## [3] "GlobBiomass" "CCI"
## [5] "This study (cap999pctl_maskWUwb)"
# Prep DF
df_lc <- stats_by_lc %>%
filter(!is.na(name)) %>%
mutate(pct_of_lc_NA = 1 - pct_of_lc_100m) %>%
pivot_longer(cols = starts_with('pct_of_lc'),
names_prefix = 'pct_of_lc_',
names_to = 'group',
values_to = 'pct') %>%
mutate(group = factor(group, levels = c('NA', '100m')),
name = factor(name, levels = map_order),
LC = recode_factor(LC, !!!lc_tbl))
# Plot
(bp <- df_lc %>%
ggplot(aes(x = pct, y = name, fill = group)) +
geom_col(position = 'fill', show.legend = FALSE) +
scale_x_continuous(n.breaks = 3) +
theme_minimal() +
labs(fill="",x=NULL,y=NULL,title="",caption="") +
facet_wrap(facets = vars(LC)))
# Save
ggsave(here::here(comparison_dir, '07_pct_masked_by_lc.png'),
width = 6, height = 4, dpi = 120)
# Get percent of total
pct_total <- df_lc %>%
group_by(name, map_code) %>%
summarize(n_vals_100m = sum(n_vals_100m, na.rm = TRUE),
n_lc_100m = sum(n_lc_100m, na.rm = TRUE)) %>%
dplyr::mutate(pct_of_lc_100m = n_vals_100m / n_lc_100m,
pct_of_lc_NA = 1 - pct_of_lc_100m) %>%
pivot_longer(cols = starts_with('pct_of_lc'),
names_prefix = 'pct_of_lc_',
names_to = 'group',
values_to = 'pct') %>%
dplyr::mutate(group = factor(group, levels = c('NA', '100m')))
# Plot
(bp <- pct_total %>%
ggplot(aes(x = pct, y = name, fill = group)) +
geom_col(position = 'fill', show.legend = FALSE) +
scale_x_continuous(n.breaks = 3) +
theme_minimal() +
labs(fill="",x=NULL,y=NULL,title="",caption="") +
ggtitle('Percent of all pixels (aggregated to 100m)'))
# Table
(d2 <- pct_total %>%
dplyr::mutate(pct_coverage = str_c(round(pct*100, digits = 1), '%')) %>%
dplyr::arrange(desc(pct)) %>%
dplyr::filter(!group %in% c('NA')) %>%
dplyr::select(name, pct_coverage))
tab <- gridExtra::tableGrob(d2, rows = NULL, cols = c('Map', 'Percent of land'))
# Display side-by-side
wrap_plots(tab, bp)
# Scatterplot - AGB against backscatter ----------------------------------------
plot_scatter_ext <- function(plot_means, y_var, name){
df <- plot_means %>% transmute(diff = .data[[y_var]] - AGB_ha) %>% deframe
ggplot(plot_means, aes(x=AGB_ha, y=.data[[y_var]])) +
geom_point() +
labs(x = expression(paste("Observed AGB (Mg ha"^"-1", ")")),
y = expression(paste("Estimated AGB (Mg ha"^"-1", ")")),
subtitle = str_c(name, ", mean difference: ",
round(mean(df, na.rm=T), 1), " +/- ",
round(sd(df, na.rm=T), 1), " t/ha")) +
coord_cartesian(xlim=c(0,160), ylim=c(0,160)) +
geom_abline(intercept = 0, slope = 1, linetype='dashed', alpha=.5) +
geom_text(label=rownames(plot_means), hjust=0, vjust=0, size=3.5) +
theme_minimal()
}
# Load data
plot_means <- read_csv(plot_ext_csv)
# Extract just AGB and backscatter as dataframe
p0 <- plot_scatter_ext(plot_means, 'internal', name='Our AGB')
p1 <- plot_scatter_ext(plot_means, 'glob', name='GlobBiomass')
p2 <- plot_scatter_ext(plot_means, 'esa', name='ESA')
p3 <- plot_scatter_ext(plot_means, 'avit', name='Avitabile')
p4 <- plot_scatter_ext(plot_means, 'bacc', name='Baccini')
# Save
plot_fp <- here::here(comparison_dir,
str_c('comparison_scatter_', agb_code, '_ourAGB.png'))
ggsave(plot_fp, p0, width=10, height=10, units='cm')
(p1 | p2 )/ (p3 | p4)
plot_fp <- here::here(comparison_dir,
str_c('comparison_scatter_', agb_code, '.png'))
ggsave(plot_fp, width=20, height=20, units='cm')
# Load
sum_tbl <- read_csv(compare_by_lc_csv) %>%
filter(!is.na(R2),
mask != 'WUBGS')
# Experiment with views
long_tbl <- sum_tbl %>%
rename(product = name) %>%
pivot_longer(R2:med_diff, names_repair = 'unique')
long_tbl %>%
filter(name %in% c('R2', 'RMSD', 'bias')) %>%
ggplot(aes(x = product, y = value, color = product)) +
geom_point() +
facet_grid(rows = vars(name), cols = vars(mask), scales = 'free_y') +
theme(axis.title = element_blank())
long_tbl %>%
filter(name %in% c('MBD')) %>%
ggplot(aes(x = product, y = value, color = product)) +
geom_point() +
geom_hline(yintercept = 0) +
facet_grid(rows = vars(name), cols = vars(mask), scales = 'free_y') +
theme(axis.title = element_blank())
long_tbl %>%
filter(name %in% c('R2', 'RMSD', 'bias'),
mask %in% c('TC', 'BGS')) %>%
ggplot(aes(x = mask, y = value, color = product)) +
geom_point() +
facet_wrap(vars(name), scales = 'free_y') +
theme(axis.title = element_blank())
# Tree cover
mskvals <- mskvals_list[[1]]
agb_fps_sub <- agb_fps[c(2,3,4,5)]
plots_tc <- agb_fps_sub %>%
purrr::map(compare_by_given_lc, agb_fp, lc_fp,
mskvals = mskvals, comparison_dir,
return_obj = 'plot',
boundary_fp = hti_poly_fp) %>%
flatten()
plots_tc %>% saveRDS(here::here(comparison_dir, 'plot_tc.rds'))
# Plots for given mask
(patches <- plots_tc %>%
wrap_plots(nrow = 4, guides = 'collect') &
theme(axis.title = element_blank()))
map_fp <- here::here(comparison_dir, str_glue("scattdens_{mskvals$code}.png"))
ggsave(map_fp, width=4, height=12, dpi=120)
# Total
mskvals <- mskvals_list[[2]]
agb_fps_sub <- agb_fps[c(2,3,4,5)]
plots_total <- agb_fps_sub %>%
purrr::map(compare_by_given_lc, agb_fp, lc_fp,
mskvals = mskvals, comparison_dir,
return_obj = 'plot',
boundary_fp = hti_poly_fp) %>%
flatten()
plots_total %>% saveRDS(here::here(comparison_dir, 'plots_total.rds'))
# Plots for given mask
(patches <- plots_total %>%
wrap_plots(nrow = 4, guides = 'collect') &
theme(axis.title = element_blank()))
map_fp <- here::here(comparison_dir, str_glue("scattdens_{mskvals$code}.png"))
ggsave(map_fp, width=4, height=12, dpi=120)
# BGS
mskvals <- mskvals_list[[3]]
agb_fps_sub <- agb_fps[c(2,3,4,5)]
plots_bgs <- agb_fps_sub %>%
purrr::map(compare_by_given_lc, agb_fp, lc_fp,
mskvals = mskvals, comparison_dir,
return_obj = 'plot',
boundary_fp = hti_poly_fp) %>%
flatten()
plots_bgs %>% saveRDS(here::here(comparison_dir, 'plots_bgs.rds'))
# Plots for given mask
(patches <- plots_bgs %>%
wrap_plots(nrow = 4) &
theme(axis.title = element_blank()))
map_fp <- here::here(comparison_dir, str_glue("scattdens_{mskvals$code}.png"))
ggsave(map_fp, width=4, height=12, dpi=120)
# Combine
(patches <- c(plots_tc, plots_bgs, plots_total) %>%
wrap_plots(nrow = 4, ncol = 3,
byrow = FALSE,
guides = 'collect') &
theme(axis.title = element_blank()))
map_fp <- here::here(comparison_dir, str_glue("scattdens_all.png"))
ggsave(map_fp, width=12, height=12, dpi=120)
plots_tc$glob +
guides(color = guide_legend('density'),
alpha = guide_legend('density'))
plots_tc$glob + guides(alpha = FALSE)
# Tree cover
mskvals <- mskvals_list[[1]]
agb_fps_sub <- agb_fps[c(2,3,4,5)]
df <- compare_by_given_lc(agb_fps_sub[[1]], agb_fp, lc_fp,
mskvals = mskvals, comparison_dir,
return_obj = 'diffs',
boundary_fp = hti_poly_fp)
ext_name = ''
sample_size = 50000
xlim = c(0, 300)
# Guides are integrated where possible
p + guides(colour = guide_legend("title"), size = guide_legend("title"),
shape = guide_legend("title"))
# ggplot object
dat <- data.frame(x = 1:5, y = 1:5, p = 1:5, q = factor(1:5),
r = factor(1:5))
p <- ggplot(dat, aes(x, y, colour = p, size = q, shape = r)) + geom_point()
# without guide specification
p
# Show colorbar guide for colour.
# All these examples below have a same effect.
p + guides(colour = "colorbar", size = "legend", shape = "legend")
p + guides(colour = guide_colorbar(), size = guide_legend(),
shape = guide_legend())
p +
scale_colour_continuous(guide = "colorbar") +
scale_size_discrete(guide = "legend") +
scale_shape(guide = "legend")
# Remove some guides
p + guides(colour = "none")
p + guides(colour = "colorbar",size = "none")
# Guides are integrated where possible
p + guides(colour = guide_legend("title"), size = guide_legend("title"),
shape = guide_legend("title"))
# same as
g <- guide_legend("title")
p + guides(colour = g, size = g, shape = g)
p + theme(legend.position = "bottom")
# position of guides
# Set order for multiple guides
ggplot(mpg, aes(displ, cty)) +
geom_point(aes(size = hwy, colour = cyl, shape = drv)) +
guides(
colour = guide_colourbar(order = 1),
shape = guide_legend(order = 2),
size = guide_legend(order = 3)
)
plot_agb_map <- function(x, mask_poly = NULL, borders = NULL, resolution = 200) {
lyr_name <- x$name
# Load AGB
agb <- rast(x$fp)
# Conditionally aggregate map based on resolution
resx <- res(agb) %>% nth(1) #%>% signif(2)
fact <- resolution/resx/113000
agg_factor <- ifelse(fact >= 1.5, round(fact, 0), 1)
if (agg_factor > 1) {
agb <- terra::aggregate(agb, fact = agg_factor, fun = 'mean', na.rm = TRUE)
} else if(fact < 1) {
print('WARNING: input map resolution is greater than target')
}
# Convert raster to dataframe
agb_a_df <- agb %>% as.data.frame(xy = T) %>% rename(agb = 3)
if(is.null(mask_poly)){
# Urban polygons
tol <- 0.001
urb_poly_fp <- here::here(tidy_dir, 'landcover', str_glue('lc2_urban_tol{tol}.gpkg'))
if(!file.exists(urb_poly_fp)) {
mask_poly <- st_read(lc_pols_fp) %>%
filter(LC == 2) %>%
st_simplify(dTolerance = tol) %>%
st_union()
mask_poly %>% st_write(urb_poly_fp)
} else {
mask_poly <- st_read(urb_poly_fp)
}
}
if(is.null(mask_poly)){
# Haiti administrative boundaries
borders <- st_read(here::here(tidy_dir, 'contextual_data', 'HTI_adm', 'HTI_adm1.shp'))
}
# Plot
agb_map <- ggplot(agb_a_df) +
geom_raster(aes(x = x, y = y, fill = agb)) +
# geom_tile(aes(x = x, y = y, fill = agb)) +
geom_sf(data = mask_poly, fill = "gray63", color = NA) +
geom_sf(data = borders, fill = NA, color = "gray39", size = 0.3) +
scale_fill_gradientn(expression(paste("Aboveground\nbiomass (Mg ha"^"-1", ")")),
na.value = "transparent",
colors = agb1_palette,
limits = c(0, 120),
oob = scales::squish
) +
ggthemes::theme_map() +
ggtitle(lyr_name) +
theme(legend.position=c(0.23, 0.9),
legend.title.align=0,
legend.justification = c(1,1),
legend.box.background = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_blank(),
plot.title = element_text(vjust = -8, hjust = 0.048))
}
# Urban polygons
tol <- 0.001
urb_poly_fp <- here::here(tidy_dir, 'landcover', str_glue('lc2_urban_tol{tol}.gpkg'))
if(!file.exists(urb_poly_fp)) {
urban <- st_read(lc_pols_fp) %>%
filter(LC == 2) %>%
st_simplify(dTolerance = tol) %>%
st_union()
urban %>% st_write(urb_poly_fp)
} else {
urban <- st_read(urb_poly_fp)
}
## Reading layer `lc2_urban_tol0.001' from data source `/home/esturdivant/PROJECTS/biomass-espanola/data/tidy/landcover/lc2_urban_tol0.001.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 0 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.4575 ymin: 18.04225 xmax: -71.68825 ymax: 20.04025
## Geodetic CRS: WGS 84
# Haiti administrative boundaries
hti_adm1 <- st_read(here::here(tidy_dir, 'contextual_data', 'HTI_adm', 'HTI_adm1.shp'))
## Reading layer `HTI_adm1' from data source `/home/esturdivant/PROJECTS/biomass-espanola/data/tidy/contextual_data/HTI_adm/HTI_adm1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 10 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.48125 ymin: 18.0218 xmax: -71.61815 ymax: 20.09042
## Geodetic CRS: WGS 84
resolution <- 400
(agb_map <- plot_agb_map(agb_fps$internal, urban, hti_adm1, resolution) +
theme(plot.title = element_blank()))
Data resampled to 400 m resolution for plotting.
# Extract plot means for all AGB maps
agb_fps$internal$name <- 'This study'
maps <- agb_fps %>%
purrr::map(plot_agb_map, urban, hti_adm1, resolution = resolution)
## [1] "WARNING: input map resolution is greater than target"
# Set output dir for side-by-side maps
out_dir <- here::here(comparison_dir, 'sidebyside')
dir.create(out_dir)
# Ours, with GlobBiomass, ESA, and Baccini
all_maps <- maps[c(1,3,2,5)] %>%
wrap_plots(nrow = 2, guides = 'collect')
map_fp <- here::here(out_dir, str_glue("sidebyside_Us_ESA_GlobB_Bacc_{resolution}m.png"))
ggsave(map_fp, width=12, height=10, dpi=120)
all_maps
# 4 external maps
all_maps <- maps[c(4,2,3,5)] %>%
wrap_plots(nrow = 2, guides = 'collect')
map_fp <- here::here(out_dir, str_glue("sidebyside_Avit_GlobB_ESA_Bacc_{resolution}m.png"))
ggsave(map_fp, width=12, height=10, dpi=120)
all_maps
# All 5, equal sizes
maps$space <- guide_area()
all_maps <- maps[c(1, 3, 2, 6, 4, 5)] %>%
wrap_plots(guides = 'collect')
map_fp <- here::here(out_dir, str_glue("all5_equal_{resolution}m.png"))
ggsave(map_fp, width=14, height=9, dpi=120)
all_maps
# All 5, with ours bigger
four <- ((maps$bacc + maps$avit) / (maps$glob + maps$esa)) &
theme(legend.position = 'none')
all_maps <- (maps$internal + plot_layout(guides = 'keep')) / four
map_fp <- here::here(out_dir, str_glue("all5_{resolution}m.png"))
ggsave(map_fp, all_maps, width=5, height=8, dpi=120)
all_maps
# ~ Difference maps (only for all LCs) -----------------------------------------
mskvals <- list(values = c(1,2,3,4,5,6), code = 'Total')
diffmaps <- agb_fps[c(2,3,4,5)] %>%
purrr::map(compare_by_given_lc, agb_fp, lc_fp,
mskvals = mskvals, comparison_dir,
return_obj = 'map',
boundary_fp = hti_poly_fp)
# 4 external maps
all_maps <- diffmaps %>%
wrap_plots(nrow = 2, guides = 'collect') &
ggthemes::theme_map()
map_fp <- here::here(out_dir, str_glue("difference_maps.png"))
ggsave(map_fp, width=14, height=10, dpi=120)
all_maps